require("Hmisc")
require("plyr")
require("ggplot2")
require("car")
require("compositions")
options(digits=10)
options(scipen=10)
#setwd("/nfs/projects/k/kyu")
setwd("U:/gov2001/paper")

##############################################################
## FIGURE 1
##############################################################

load("results_volat_all.RData")
load("results_beta_all.RData")
#load("results_volat_top1000.RData")
#load("results_beta_top1000.RData")

# Volatility plot
d <- volat
d$YEAR <- 1963+volat[volat$VOLATQ==1,"PERIOD"]/12
d$"Volatility Quintiles" <- recode(d$VOLATQ,"1='Quintile 1 (Lowest)';2='Quintile 2';3='Quintile 3';4='Quintile 4';5='Quintile 5 (Highest)'")
g <- ggplot(data=d,aes(x=YEAR,y=CUMRET,colour=`Volatility Quintiles`,linetype=`Volatility Quintiles`)) + geom_line(lwd=1)
g <- g + scale_y_log10("Cumulative Return",limits = c(.1,100),breaks=c(.1,1,10,100),labels=c("$0.10","$1","$10","$100"))
g <- g + scale_x_continuous("",limits = c(1968,2009),breaks=seq(1968,2009,5),seq(1968,2009,5),expand=c(0,0))
g <- g + ggtitle("A. All Stocks, Volatility Quintiles\nCumulative Return of $1 invested in 1968\n")
#g <- g + ggtitle("B. Top 1,000 Stocks, Volatility Quintiles\nCumulative Return of $1 invested in 1968\n")
g <- g + scale_color_manual(values=rainbow(5))
g <- g + scale_linetype_manual(values=c(1,3,3,3,1))
g <- g + theme(legend.position=c(.2,.8))
g
dev.copy(pdf,"plotA.pdf")
#dev.copy(pdf,"plotB.pdf")
dev.off()

# Beta plot
d <- beta
d$YEAR <- 1963+volat[beta$BETAQ==1,"PERIOD"]/12
d$"Beta Quintiles" <- recode(d$BETAQ,"1='Quintile 1 (Lowest)';2='Quintile 2';3='Quintile 3';4='Quintile 4';5='Quintile 5 (Highest)'")
g <- ggplot(data=d,aes(x=YEAR,y=CUMRET,colour=`Beta Quintiles`,linetype=`Beta Quintiles`)) + geom_line(lwd=1)
g <- g + scale_y_log10("Cumulative Return",limits = c(.1,100),breaks=c(.1,1,10,100),labels=c("$0.10","$1","$10","$100"))
g <- g + scale_x_continuous("",limits = c(1968,2009),breaks=seq(1968,2009,5),seq(1968,2009,5),expand=c(0,0))
g <- g + ggtitle("C. All Stocks, Beta Quintiles\nCumulative Return of $1 invested in 1968\n")
#g <- g + ggtitle("D. Top 1,000 Stocks, Beta Quintiles\nCumulative Return of $1 invested in 1968\n")
g <- g + scale_color_manual(values=rainbow(5))
g <- g + scale_linetype_manual(values=c(1,3,3,3,1))
g <- g + theme(legend.position=c(.2,.8))
g
dev.copy(pdf,"plotC.pdf")
#dev.copy(pdf,"plotD.pdf")
dev.off()

##############################################################
## TABLE 1
##############################################################

# Load data and calculate monthly returns
load("results_data_all.RData")
v <- ddply(results,.(PERIOD,VOLATQ),summarise,RETURN=weighted.mean(RETURN,CAP))
b <- ddply(results,.(PERIOD,BETAQ),summarise,RETURN=weighted.mean(RETURN,CAP))
load("results_data_top1000.RData")
v1000 <- ddply(results,.(PERIOD,VOLATQ),summarise,RETURN=weighted.mean(RETURN,CAP))
b1000 <- ddply(results,.(PERIOD,BETAQ),summarise,RETURN=weighted.mean(RETURN,CAP))
data.rf <- read.delim("famafrench.csv",colClasses=c("numeric","numeric","numeric","numeric","numeric"))
rf <- data.rf[data.rf$DATE>=196801,"RF"]
rmrf <- data.rf[data.rf$DATE>=196801,"RMRF"]

# Calculate cumulative returns
start <- 61
end <- 552
cumret.v <- ddply(v,.(VOLATQ),summarise,CUMRET=prod(RETURN+1))[,2]
cumret.v <- c(cumret.v,ddply(v1000,.(VOLATQ),summarise,CUMRET=prod(RETURN+1))[,2])
cumret.b <- ddply(b,.(BETAQ),summarise,CUMRET=prod(RETURN+1))[,2]
cumret.b <- c(cumret.b,ddply(b1000,.(BETAQ),summarise,CUMRET=prod(RETURN+1))[,2])

#################
# Volatilty sorts
#################

# Arrange all stock quintiles as columns 1-5
t.out <- NA
for(i in 1:5) t.out <- cbind(t.out,v[v$VOLATQ==i,"RETURN"])
t.out <- t.out[,-1]
# Arrange top 1,000 stock quintiles as columns 6-10 
for(i in 1:5) t.out <- cbind(t.out,v1000[v1000$VOLATQ==i,"RETURN"])
# Multiply by 100
t.out <- t.out*100
# Minus risk-free rate or market return
t.rf <- t.out-rf
rm <- rmrf+rf
t.rm <- t.out-rm

# Geometric Mean (rf)
geom <- (apply(100+t.rf,2,geometricmean)-100)*12
# Arithmetic Mean (rf)
meanrf <- apply(t.rf,2,mean)*12
# SD
stdev <- apply(t.out,2,sd)*sqrt(12)
# Sharpe ratio
sharpe <- meanrf/stdev
# Arithmetic Mean (rm)
meanrm <- apply(t.rm,2,mean)*12
# Tracking error
te <- sqrt(apply((t.out-(rm-rf))^2,2,mean))*sqrt(12)
# Information ratio
ir <- meanrm/te
# Beta
pbeta <- apply(t.out,2,function(x){coef(lm(x~rmrf))[2]})
# Alpha
mat.rf <- matrix(rep(rf,10),nrow=492,ncol=10)
mat.rmrf <- t(t(matrix(rep(rmrf,10),nrow=492,ncol=10))*pbeta)
mat.alpha <- t.out-(mat.rf+mat.rmrf)
alpha <- apply(mat.alpha,2,mean)*12

# Assemble table
table1a <- rbind(cumret.v,geom,meanrf,stdev,sharpe,meanrm,te,ir,pbeta,alpha)

#################
# Beta sorts
#################

# Arrange all stock quintiles as columns 1-5
t.out <- NA
for(i in 1:5) t.out <- cbind(t.out,b[b$BETAQ==i,"RETURN"])
t.out <- t.out[,-1]
# Arrange top 1,000 stock quintiles as columns 6-10 
for(i in 1:5) t.out <- cbind(t.out,b1000[b1000$BETAQ==i,"RETURN"])
# Multiply by 100
t.out <- t.out*100
# Minus risk-free rate or market return
t.rf <- t.out-rf
rm <- rmrf+rf
t.rm <- t.out-rm

# Geometric Mean (rf)
geom <- (apply(100+t.rf,2,geometricmean)-100)*12
# Arithmetic Mean (rf)
meanrf <- apply(t.rf,2,mean)*12
# SD
stdev <- apply(t.out,2,sd)*sqrt(12)
# Sharpe ratio
sharpe <- meanrf/stdev
# Arithmetic Mean (rm)
meanrm <- apply(t.rm,2,mean)*12
# Tracking error
te <- sqrt(apply((t.out-(rm-rf))^2,2,mean))*sqrt(12)
# Information ratio
ir <- meanrm/te
# Beta
pbeta <- apply(t.out,2,function(x){coef(lm(x~rmrf))[2]})
# Alpha
mat.rf <- matrix(rep(rf,10),nrow=492,ncol=10)
mat.rmrf <- t(t(matrix(rep(rmrf,10),nrow=492,ncol=10))*pbeta)
mat.alpha <- t.out-(mat.rf+mat.rmrf)
alpha <- apply(mat.alpha,2,mean)*12

# Assemble table
table1b <- rbind(cumret.b,geom,meanrf,stdev,sharpe,meanrm,te,ir,pbeta,alpha)

# Output
table1a
table1b











